#How Do Observers Assess Resolve?
#July 7, 2018
#BJPS Replication 8.R

#This file contains the replication code needed to replicate the analysis from Supplementary Appendix §1.5, 4, and 5 (response latency results, sample and weighting information, and elite experimental benchmarks)

library(stargazer)
library(ggplot2)

### To generate the necessary dataframes, either run BJPS Replication 1.R, or load the following saved objects:
library(here) #Avoids needing to setwd()
dat3 <- get(load(file="Resolve Conjoint Cleaned No US 070718.RData"))
dat_id <- get(load(file="Resolve Demographics 070718.RData"))
dat.elite <- get(load(file="Elite Benchmarks.RData"))

#Load parallel processing and recode functions in case you haven't run BJPS Replication 2.R - if not, you can skip to line 28
library(snow)
cl <- makeCluster(8,"SOCK")
clusterBootS2 <- function(data){
	i <- sample(unique(data$ID_numeric),length(unique(data$ID_numeric)),replace=TRUE)
	row.nums <- NULL
	for (j in 1:length(i)){
		row.nums <- c(row.nums, which(data$ID_numeric==i[j]))
	}
	return(data[row.nums,])
}

###### Response latency effects

set.seed(43215)
B <- 1500

#First, create new dataset with measures of whether each treatment combination in a given round is conjunctive vs disjunctive (e.g. two profiles with the same level of capabilities, rather than one with high capabilities and the other with low)

datR <- data.frame(expand.grid(ID_numeric=unique(dat3$ID_numeric),Round=c(1:8), responseTime=NA, dRegimeType=NA, dCapabilities=NA, dStakes=NA, dNewLeader=NA, dMilService=NA, dmaleLeader=NA, dactor=NA, dcurBehavior=NA, dprevOpponent=NA, dprevRole=NA, dprevAct=NA, dprevLdr=NA))

for (i in 1:nrow(datR)){
	loc <- which(dat3$ID_numeric==datR$ID_numeric[i] & dat3$Round==datR$Round[i]) #Find in original dataset
	datR$responseTime[i]<-dat3$responseTime[loc[1]]
	temp <- dat3[loc[1],c(5:11,13,17:20)]==dat3[loc[2],c(5:11,13,17:20)] #Compare treatments across each country in a given round
	datR[i,4:15] <- as.numeric(temp)		
}
#Note that 0 values = disjunctive, 1 = conjunctive

#Create logged version of response time variable
datR$responseTime <- log(as.numeric(datR$responseTime)) 

#Now, bootstrap regression model
bootConjoint <- function(...){
	temp <- clusterBootS2(datR)
	mod.temp <- lm(responseTime ~  dCapabilities + dStakes + dRegimeType + dactor + dNewLeader + dMilService + dmaleLeader + dprevRole + dprevOpponent + dprevAct + dprevLdr + dcurBehavior, data=temp)
	return(coef(mod.temp))
}

clusterExport(cl, list("datR", "bootConjoint", "clusterBootS2"))

system.time(boot.responseTimes <- parSapply(cl, rep(1,1500), bootConjoint))

#Extract quantities of interest
RT.mat <- cbind(apply(boot.responseTimes, 1, mean), apply(boot.responseTimes, 1, quantile, c(0.025)), apply(boot.responseTimes, 1, quantile, c(0.05)), apply(boot.responseTimes, 1, quantile, c(0.95)), apply(boot.responseTimes,1,quantile, c(0.975)))

#Label factors
attributes <- c("Capabilities", "Stakes", "Regime Type", "Actor", "New Leader", "Mil Service", "Male",  "Prev Role", "Prev Opponent", "Prev Act", "Prev Leader", "Cur Behavior")

#### Figure 14: Response latency effects ####
dev.new(height=5, width=4.8)
par(oma=c(0,0,0,0), mar=c(4,0,1,0))
plot(RT.mat[,1], nrow(RT.mat):1, type="n", axes=FALSE, xlab="Effect of conjunctive cues on logged response time", ylab="",xlim=c(-0.2,0.1), ylim=c(1,12))

for (i in 1:nrow(RT.mat)){
	points(RT.mat[i,1], nrow(RT.mat)-(i-1), pch=16)
	lines(RT.mat[i,c(2,5)], rep(nrow(RT.mat)-(i-1),2))
	lines(RT.mat[i,c(3,4)], rep(nrow(RT.mat)-(i-1),2), lwd=2)		
	text(-0.15, nrow(RT.mat)-(i), attributes[i], cex=0.7)			
	}
abline(v=0)
axis(1)
box()

###### Sample characteristics

#Let's set the weight at 0 for the one missing value
dat_id$eWeight[which(is.na(dat_id$eWeight))] <- 0

tab4 <- cbind(c(round(mean(dat_id$male, na.rm=TRUE), digits=3), round(mean(dat_id$age1, na.rm=TRUE), digits=3), round(mean(dat_id$age2, na.rm=TRUE), digits=3), round(mean(dat_id$age3, na.rm=TRUE), digits=3), round(mean(dat_id$age4, na.rm=TRUE), digits=3), round(mean(dat_id$educ1, na.rm=TRUE), digits=3), round(mean(dat_id$educ2, na.rm=TRUE), digits=3), round(mean(dat_id$educ3, na.rm=TRUE), digits=3), round(mean(dat_id$educ4, na.rm=TRUE), digits=3)), c(round(weighted.mean(dat_id$male, dat_id$eWeight, na.rm=TRUE), digits=3), round(weighted.mean(dat_id$age1, dat_id$eWeight, na.rm=TRUE), digits=3), round(weighted.mean(dat_id$age2, dat_id$eWeight, na.rm=TRUE), digits=3), round(weighted.mean(dat_id$age3, dat_id$eWeight, na.rm=TRUE), digits=3), round(weighted.mean(dat_id$age4, dat_id$eWeight, na.rm=TRUE), digits=3), round(weighted.mean(dat_id$educ1, dat_id$eWeight, na.rm=TRUE), digits=3), round(weighted.mean(dat_id$educ2, dat_id$eWeight, na.rm=TRUE), digits=3), round(weighted.mean(dat_id$educ3, dat_id$eWeight, na.rm=TRUE), digits=3), round(weighted.mean(dat_id$educ4, dat_id$eWeight, na.rm=TRUE), digits=3)), c(0.492, 0.128, 0.342, 0.341, 0.189, 0.420, 0.194, 0.282, 0.104))

rownames(tab4) <- c("Male", "18 to 24 years", "25 to 44 years", "45 to 64 years", "65 years and over", "High School or less", "Some college", "College/University", "Graduate/Professional school")
colnames(tab4) <- c("Unweighted Sample", "Weighted Sample", "Population Target")


#### Table 4: Sample characteristics
stargazer(tab4)

######### Appendix §5: elite benchmarks

#Note: In order to protect the anonymity of our Knesset respondents, and consistent with IRB rules about removing identifiable information from publically-available data, a small number of demographic characteristics that would render respondents identifiable have been removed from the Knesset and population frame data (e.g. gender, age, number of terms in Knesset, cabinet member experience, and political party).  Any analysis that uses these variables has been flagged in the replication code below, and modified models included where applicable; for the raw R output, see Tables A5-A7 Output.pdf

##### Table 5: Knesset Sample Characteristics
#Knesset Member:
round(table(dat.elite$currentMK)/sum(table(dat.elite$currentMK)),digits=2) #25% current, 75% former
#Exp on Foreign Affairs/Defense Committee
round(table(dat.elite$FAexp1)/sum(table(dat.elite$FAexp1)),digits=2) #67% backup or full
round(table(dat.elite$FAexp2)/sum(table(dat.elite$FAexp2)),digits=2) #54% full
#Highest level of experience
#round(table(dat.elite$eliteExperience)/sum(table(dat.elite$eliteExperience)),digits=2) #58% not a minister, 29% deputy minister, 12% cabinet minister or higher
#Male
#round(table(dat.elite$male)/sum(table(dat.elite$male)),digits=2) #84%
#Served in military
#army <- na.omit(as.numeric(dat.elite$army))
#round((table(army)[2] + table(army)[3]) /sum(table(army)),digits=2) #95% military service
#round((table(army)[3]) /sum(table(army)),digits=2) #64.3% combat
#Age
#round(mean(dat.elite$age,na.rm=TRUE),digits=1) #61.4
#round(sd(dat.elite$age,na.rm=TRUE),digits=1) #10.7
#Terms in Knesset
#round(mean(dat.elite$terms,na.rm=TRUE),digits=1) #3.0
#round(sd(dat.elite$terms,na.rm=TRUE),digits=1) #2.1
#Military assertiveness
round(mean(dat.elite$milAssert1,na.rm=TRUE),digits=2) #0.61
round(sd(dat.elite$milAssert1,na.rm=TRUE),digits=2) #0.20
#Right wing ideology
round(mean(dat.elite$ideo1,na.rm=TRUE),digits=2) #0.45
round(sd(dat.elite$ideo1,na.rm=TRUE),digits=2) #0.24
#Hardline (Arab-Israeli)
round(mean(dat.elite$hawk,na.rm=TRUE),digits=2) #0.39
round(sd(dat.elite$hawk,na.rm=TRUE),digits=2) #0.25
#International trust
round(mean(dat.elite$intTrust1,na.rm=TRUE),digits=2) #0.40
round(sd(dat.elite$intTrust1,na.rm=TRUE),digits=2) #0.25

#Missingness information for footnote 7
table(is.na(dat.elite$currentMK))
table(is.na(dat.elite$FAexp1))
table(is.na(dat.elite$FAexp2))
#table(is.na(dat.elite$eliteExperience))
#table(is.na(dat.elite$male))
#table(is.na(dat.elite$army)) #5 missing
#table(is.na(dat.elite$age))
#table(is.na(dat.elite$terms))
table(is.na(dat.elite$milAssert1))
table(is.na(dat.elite$ideo1)) #2 missing
table(is.na(dat.elite$hawk)) #2 missing
table(is.na(dat.elite$intTrust1)) #3 missing
table(is.na(dat.elite$outcmBW)) #0 missing
table(is.na(dat.elite$DV_A1)) #1 missing

##### Table 6: Sample representativeness tests

#pop.elite2 <- subset(pop.elite, pop.elite$Contacted==1) #Create dataframe for contacted MKs only

#mod.rep1 <- lm(Participated ~ currentMK  + as.factor(eliteExperience) + male + terms, data=pop.elite)
#mod.rep1a <- lm(Participated ~ currentMK  + as.factor(eliteExperience) + male + terms + LRscale, data=pop.elite)
#mod.rep2 <- lm(Participated ~ currentMK  + as.factor(eliteExperience) + male + terms, data=pop.elite2)
#mod.rep2a <- lm(Participated ~ currentMK  + as.factor(eliteExperience) + male + terms + LRscale, data=pop.elite2)

#stargazer(mod.rep1, mod.rep1a, mod.rep2, mod.rep2a, omit.stat=c("LL", "ser", "f"), style="apsr", digits=3, label="tab:a6")

##### Table 7: Elite experiment balance checks

randMat <- cbind(round(c(mean(dat.elite$currentMK[which(dat.elite$A==0)], na.rm=TRUE), mean(dat.elite$FAexp1[which(dat.elite$A==0)], na.rm=TRUE), mean(dat.elite$eliteExperience[which(dat.elite$A==0)], na.rm=TRUE), mean(dat.elite$male[which(dat.elite$A==0)], na.rm=TRUE), mean(dat.elite$age[which(dat.elite$A==0)], na.rm=TRUE), mean(dat.elite$combat1[which(dat.elite$A==0)], na.rm=TRUE), mean(dat.elite$milAssert1[which(dat.elite$A==0)], na.rm=TRUE), mean(dat.elite$ideo1[which(dat.elite$A==0)], na.rm=TRUE), mean(dat.elite$hawk[which(dat.elite$A==0)], na.rm=TRUE), mean(dat.elite$intTrust1[which(dat.elite$A==0)], na.rm=TRUE)), digits=2), round(c(mean(dat.elite$currentMK[which(dat.elite$A==1)], na.rm=TRUE), mean(dat.elite$FAexp1[which(dat.elite$A==1)], na.rm=TRUE), mean(dat.elite$eliteExperience[which(dat.elite$A==1)], na.rm=TRUE), mean(dat.elite$male[which(dat.elite$A==1)], na.rm=TRUE), mean(dat.elite$age[which(dat.elite$A==1)], na.rm=TRUE), mean(dat.elite$combat1[which(dat.elite$A==1)], na.rm=TRUE), mean(dat.elite$milAssert1[which(dat.elite$A==1)], na.rm=TRUE), mean(dat.elite$ideo1[which(dat.elite$A==1)], na.rm=TRUE), mean(dat.elite$hawk[which(dat.elite$A==1)], na.rm=TRUE), mean(dat.elite$intTrust1[which(dat.elite$A==1)], na.rm=TRUE)), digits=2), round(c(mean(dat.elite$currentMK[which(dat.elite$Threat==0)], na.rm=TRUE), mean(dat.elite$FAexp1[which(dat.elite$Threat==0)], na.rm=TRUE), mean(dat.elite$eliteExperience[which(dat.elite$Threat==0)], na.rm=TRUE), mean(dat.elite$male[which(dat.elite$Threat==0)], na.rm=TRUE), mean(dat.elite$age[which(dat.elite$Threat==0)], na.rm=TRUE), mean(dat.elite$combat1[which(dat.elite$Threat==0)], na.rm=TRUE), mean(dat.elite$milAssert1[which(dat.elite$Threat==0)], na.rm=TRUE), mean(dat.elite$ideo1[which(dat.elite$Threat==0)], na.rm=TRUE), mean(dat.elite$hawk[which(dat.elite$Threat==0)], na.rm=TRUE), mean(dat.elite$intTrust1[which(dat.elite$Threat==0)], na.rm=TRUE)), digits=2), round(c(mean(dat.elite$currentMK[which(dat.elite$Threat==1)], na.rm=TRUE), mean(dat.elite$FAexp1[which(dat.elite$Threat==1)], na.rm=TRUE), mean(dat.elite$eliteExperience[which(dat.elite$Threat==1)], na.rm=TRUE), mean(dat.elite$male[which(dat.elite$Threat==1)], na.rm=TRUE), mean(dat.elite$age[which(dat.elite$Threat==1)], na.rm=TRUE), mean(dat.elite$combat1[which(dat.elite$Threat==1)], na.rm=TRUE), mean(dat.elite$milAssert1[which(dat.elite$Threat==1)], na.rm=TRUE), mean(dat.elite$ideo1[which(dat.elite$Threat==1)], na.rm=TRUE), mean(dat.elite$hawk[which(dat.elite$Threat==1)], na.rm=TRUE), mean(dat.elite$intTrust1[which(dat.elite$Threat==1)], na.rm=TRUE)), digits=2)) #Will give warnings for highest level of experience, male, age, which have been removed from the data

rownames(randMat) <- c("Current member", "Foreign affairs experience", "Highest level of experience", "Male", "Age", "Active combat experience", "Military assertiveness", "Right wing ideology", "Hawkishness (Arab-Israeli conflict)", "International Trust")

stargazer(randMat, digits=2)

###############################
###### Downsampling mass public results

#Before we downsample the mass public results, let's recalculate our elite benchmarks

#Regime type results from elite sample
set.seed(43215)
knes.dat1 <- na.omit(as.data.frame(cbind(DV_A1=dat.elite$DV_A1, A=dat.elite$A))) #Drop missingness
boot.k1 <- rep(NA, B)
for (i in 1:B){
  j <- sample(1:nrow(knes.dat1), replace=TRUE)
  temp <- knes.dat1[j,]
  boot.k1[i] <- mean(temp$DV_A1[temp$A==1], na.rm=TRUE) - mean(temp$DV_A1[temp$A==0], na.rm=TRUE) 
}

#Costly signaling results from elite sample
set.seed(43215)
knes.dat2 <- na.omit(as.data.frame(cbind(DV_B1=dat.elite$outcmBW, Threat=dat.elite$Threat, Mobilize=dat.elite$Mobilize))) #Drop missingness
boot.k2 <- matrix(NA, B, ncol=2)
for (i in 1:B){
  j <- sample(1:nrow(knes.dat2), replace=TRUE)
  temp <- knes.dat2[j,]
  boot.k2[i,1] <- mean(temp$DV_B1[temp$Threat==1], na.rm=TRUE) 
  boot.k2[i,2] <- mean(temp$DV_B1[temp$Mobilize==1], na.rm=TRUE) 
}

#Now prepare mass data for regime type results (no mixed regimes)
conj2 <- subset(dat3, dat3$RegimeType!="Mixed")

#Regime type clustered bootstrap functions
clusterBootDS <- function(data){
	i <- sample(unique(data$ID_numeric),89,replace=TRUE) #Sample only 89 observations
	row.nums <- NULL
	for (j in 1:length(i)){
		row.nums <- c(row.nums, which(data$ID_numeric==i[j]))
	}
	return(data[row.nums,])
}

bootConjoint <- function(...){
	temp <- clusterBootDS(conj2)
	temp.res <- mean(temp$Chosen[which(temp$RegimeType=="Democracy")], na.rm=TRUE) - mean(temp$Chosen[which(temp$RegimeType=="Dictatorship")],na.rm=TRUE)
	return(temp.res)
}

clusterExport(cl, list("conj2", "bootConjoint", "clusterBootDS"))

system.time(conj.resultsD <- parSapply(cl, rep(1,B), bootConjoint)) 
conj.resultsD <- conj.resultsD * 100 #Place on the same scale as the other results

#Costly signal clustered bootstrap functions
clusterBootDS <- function(data){
	i <- sample(unique(data$ID_numeric),135,replace=TRUE) #Now sample 135 observations (so it's still 45 per cell, same as in Knesset)
	row.nums <- NULL
	for (j in 1:length(i)){
		row.nums <- c(row.nums, which(data$ID_numeric==i[j]))
	}
	return(data[row.nums,])
}

bootConjoint <- function(...){
	temp <- clusterBootDS(dat3)
	temp.res <- c(mean(temp$Chosen[which(temp$curBehavior=="Public threat")], na.rm=TRUE) - mean(temp$Chosen[which(temp$curBehavior=="Nothing")],na.rm=TRUE), mean(temp$Chosen[which(temp$curBehavior=="Mobilized troops")], na.rm=TRUE) - mean(temp$Chosen[which(temp$curBehavior=="Nothing")],na.rm=TRUE))
	return(temp.res)
}

clusterExport(cl, list("dat3", "bootConjoint", "clusterBootDS"))

system.time(conj.resultsD2 <- parSapply(cl, rep(1,B), bootConjoint)) 
conj.resultsD2 <- conj.resultsD2 * 100

#Now, combine dataframes
x_D <- data.frame(y=c(boot.k1, boot.k2[,1], boot.k2[,2], conj.resultsD, conj.resultsD2[1,], conj.resultsD2[2,]), z=rep(c("Leaders", "Public"), each=(B*3)), dv=rep(c(rep(c("Regime Type"),B), rep(c("Costly Signals"),2*B)),2), trt=rep(c(rep(c("Democracy"),B), rep(c("Public Threat"),B), rep(c("Mobilization"),B)),2))

#Collapse sample and treatment into a single variable
x_D$ztrt <- NA
x_D$ztrt[which(x_D$z=="Knesset" & x_D$trt=="Democracy")] <- "Knesset: Democracy"
x_D$ztrt[which(x_D$z=="US MTurk" & x_D$trt=="Democracy")] <- "US MTurk: Democracy"
x_D$ztrt[which(x_D$z=="Knesset" & x_D$trt=="Public Threat")] <- "Knesset: Public Threat"
x_D$ztrt[which(x_D$z=="US MTurk" & x_D$trt=="Public Threat")] <- "US MTurk: Public Threat"
x_D$ztrt[which(x_D$z=="Knesset" & x_D$trt=="Mobilization")] <- "Knesset: Mobilization"
x_D$ztrt[which(x_D$z=="US MTurk" & x_D$trt=="Mobilization")] <- "US MTurk: Mobilization"

##### Figure 15: Comparing our downsampled mass public results 
dev.new(height=3.5,width=10)
ggplot(x_D, aes(x=y)) + geom_density(aes(fill=z), alpha=0.75) + labs(x="Average Treatment Effect", y="Density") + scale_fill_brewer(guide=guide_legend(title="Sample"), palette="Set1") + facet_grid(~trt, scales="fixed") + theme_bw() + scale_y_continuous(breaks=NULL) + geom_vline(xintercept=0, linetype=3, alpha=0.5, show.legend=FALSE) 